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ABSTRACT 



Six vertical discretization schemes are compared for a 
linear, baroclinic, vorticity-divergence equation model. 
Variables are defined for two staggered and one unstaggered 
grids. -A finite difference and a Galerkin finite element 
approximation are formulated for each of the three grids. 

The models are tested in three experiments. The largest 
difference between the grids occurs in the mid-atmosphere 
diabatic heat source experiment and the unstaggered grid 
provides the best results. For the staggered grids, the 
results of the finite element models are not more accurate 
than the corresponding finite difference results. Oscilla- 
tions occur in the temperature profiles near the surface for 
both staggered finite element models. 
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I. 



INTRODUCTION 



The finite elements method (FEM) has been used in several 
atmospheric prediction models (Staniforth and Mitchell, 1977; 
Williams and Schoenstadt, 1980, Beland, et al., 1983). The 
FEM is a special case of the Galerkin procedure in which the 
dependent variables are approximated by a finite sum of 
spatially varying basis functions with time dependent coeffi- 
cients. The FEM basis functions are low order polynomials 
which are zero except in a localized region. The Galerkin 
procedure produces a set of coupled ordinary differential 
equations which are solved by using finite differences in 
time (Haltiner and Williams, 1980). 

FEM models are potentially more accurate than finite 
difference method (FDM) models. Winninghoff (1968), Arakawa 
and Lamb (1977) , and Schoenstadt (1980) have demonstrated 
the advantages of spatial staggering of predictive variables 
in finite difference approximations of shallow water equa- 
tions. The results of Williams and Schoenstadt (1980) indi- 
cated that FEM models should either use staggered nodal points 
in the momentum equations or unstaggered nodal points in the 
vorticity-divergence equations. 

The purpose of this research is to compare six linear, 
baroclinic, vorticity-divergence equation models. A finite 
difference and a finite element model are written for each of 
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the three vertical grids depicted in Figure 1. Grid A, which 
was originally developed by Lorenz (1960) for energy conser- 
vation, is a widely used grid for finite difference models. 
However Tokioka (1978) has shown that this grid has a 
computational mode in the temperature field, and Arakawa (1984) 
has found a false small scale baroclinic instability with this 
grid. Tokioka (1978) also found that grid C has computational 
modes in all fields. Beland, et al., (1983), analyzed the 
finite element formulation for grid C and found noise gener- 
ated by certain forms of friction. Grid B was introduced by 
Charney and Phillips (1953) for a finite difference quasi- 
geostrophic model. Tokioka 's (1978) analysis has shown that 
this grid has no computational modes. It is hoped that the 
finite element formulation for this grid will have high 
accuracy without the problems the other grids have had with 
computational modes. 

Each of the models used in this thesis are derived from 
the governing equations described in Chapter II. The models 
are one-wave spectral in the horizontal and have a fixed lid 
at the top of the atmosphere. Diabatic heating and mountain 
topography can be included in each model experiment. 

The results of three experiments with analytic initial 
conditions will be examined. As described in Chapter III, the 
experiments are 1) an initial perturbation in the meridional 
flow, 2) flow over mountains, and 3) flow with a mid-atmosphere 
diabatic heat source. Each experiment will be presented with 
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Figure 1. Three vertical grids to be compared for N- 
layer models. Perturbation variables are 
defined at either unstaggered levels (solid 
lines at height Z') or staggered levels 
(dashed lines at height Z) . ? is vorticity, 

D is divergence, T is potential temperature, 

T is basic state potential temperature, 

<J> is geopotential, <j> ? is surface geopotential, 
u is east-west velocity, u is basic state 
east-west velocity, v is north-south velocity, 
Z is vertical velocity, Q is dibatic heating, 
and MTS is forced vertical velocity. 
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six and sixty vertical layers. The results of the 60-layer 
models will show if the models are converging to the same 
solution. For each model, comparison of the six and 60- 
layer results indicate how well the lower resolution model 
approximates the convergent solution. The purpose of these 
experiments is to identify any problem areas in the models 
and to demonstrate the general characteristics of each model. 
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II. MODEL DESCRIPTIONS 



A. MODEL FEATURES 

The six numerical models are developed with several fea- 
tures which make the models easy to modify for a wide range 
of experiments. The depth of the atmosphere, number of 
levels and the depth of each layer are variable in each model. 
Diabatic heating and forced vertical velocity due to mountain 
topography are included in the governing equations. The user 
can prescribe heating or forced vertical velocity functions, 
or run the model with these terms defined as zero. The models 
are written in modular structure using FORTRAN '77. There is 
parallel construction between models. The subroutines used 
in one model are very similar to those used in the other five 
models. The models run very quickly on an IBM-3033 mainframe. 
A 96-hour forecast for a 12-layer FEM uses less than five 
seconds of computer processing time. 

B. GOVERNING EQUATIONS 

Each model approximates the same set of governing equa- 
tions. The vorticity (2-1), divergence (2-2), surface geo- 
potential (2-3) equations and the first law of thermodynamics 
(2-4) are the prognostic equations for the forecast variables 
vorticity, divergence, surface geopotential and potential 
temperature. The surface geopotential equation is the lower 
boundary condition on the vertical velocity. The vertical 
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coordinate Z = -ln(p/p Q ) is used, but the non-Boussinesq 



0 is the geopotential, <p = gZ, 

<P s is the surface geopotential, 

T is the potential temperature, 
u is the x-component of velocity, 
v is the y-component of velocity, 

Z is the vertical velocity, 

Q is the diabatic heating per unit time per unit mass, 

MTS is the forced vertical velocity due to flow over 
mountain topography, discussed in Chapter III, 



_ ^ 

terms involving e are replaced by one. The prognostic 
equations in the coordinates x,y,Z,t are 




0 



( 2 - 1 ) 



do . 3u. 2 dv 2 + + 

dt l 9x ; ' 3y 9y 9x 9x 3Z 




0 



( 2 - 2 ) 



d* s 



MTS 



(2-3) 



dt 



dT 

dt 



Q 



(2-4) 



Here 
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f is the Coriolis parameter, 

q . df 

S 1S 37- 

d d 3 3*3 

It = It + U 3^ + V 37 + 2 3Z' 

l ' ' 

V is the horizontal Laplacian operator. 



The prognostic equations are linearized by expanding the 
variables into the following mean state and perturbations: 



u(x,Z,t) = u(Z) + u'(x,Z,t) , 



( 2 - 5 ) 



v (x , Z , t) = v' (x , Z , t) , 



( 2 - 6 ) 



Z (x , Z , t) = Z (x , Z , t ) , 



( 2 - 7 ) 



T (x , y , Z , t ) = T(y,Z) +T'(x,Z,t) , 



( 2 - 8 ) 



(x,y, Z, t) = <My,Z) + <j>'(x,Z,t) , 



( 2 - 9 ) 



4> s (x,y , t) 



4> S (Y) + 4>g(X/t) 



( 2 - 10 ) 



C (x , Z , t) = c' (x , Z , t) , 



( 2 - 11 ) 



D (x, Z , t) = D' (x, Z, t) , 



( 2 - 12 ) 



Q (x, Z , t) = Q' (x , Z , t) , 



( 2 - 13 ) 



MTS (x , t) = MTS ' (x , t ) . 



( 2 - 14 ) 
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3v ' 



The vorticity and divergence in this system are 5 = -r — , 

0 X 

9 vi * 

D = -r — . The linearized forecast equations are 

o X 



|il = . - fD . - u f£- - 8 V , 



(2-is; 



3D' 

3 1 



f ,. . « |21- eu 

^ 3x 



, 3Z du _ d 2 <t>' 

3x dZ „ 2 ' 

3x 



(2-16) 



3 T ' 


- 3 T' , 


3 T 


; 3 t 




3 1 


- u ~ - V 

3 x 


3 y 


Z 3 Z 


+ Q' , 


34 ) ’ 


_ 


34 > 


# 




O 

3 1 


- U 3 x - V ’ 


3 y 


- RTZ 


+ MTS' , 



(2-17) 



(2-18; 



where R is the gas constant for air. 

The diagnostic variables, u', v' , Z', <p ’ , are calculated 
from the forecast variables using the definitions of divergence, 
vorticity, the hydrostatic equation and the continuity equa- 
tion. The relations are 



3u 1 
3x 



D' 



(2-19) 



3v ' 
3x 



V 



( 2 - 20 ) 



30' 

3x 



RT' 



( 2 - 21 ) 



D' 




0 . 



( 2 - 22 ) 



The use of primes to denote perturbation quantities will be 
discontinued. All quantities used in the remainder of 
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the paper will be perturbation quantities unless otherwise 
noted. 

The mean state is assumed to be in hydrostatic and geo- 
strophic balance. The term 9T/9y in the first law of thermo- 
dynamics can be evaluated by taking 3/3y of the hydrostatic 
equation and substituting for d<p/dy from the geostrophic 
relation, 3<j>/3y = -fu. Thus, 



9T f 3u 

dy R 3Z * 



(2-23) 



Geostrophic balance of the mean state at the surface implies 



= - fu S fc • 



(2-24) 



The expressions (2-23) and (2-24) are substituted into equa- 
tions (2-17) and (2-18) , respectively. 

A singlewave spectral representation is used in the 
x-direction, with wave number y = 2 it/L, where L is the wave- 
length in the x-direction. The perturbation quantities have 
the form 



£'(x,Z,t) = A^(Z,t)cos yx + A 2 (Z,t)sin yx , (2-25) 

D'(x,Z,t) = D^(Z,t)cos yx + D 2 (Z,t)sin yx , (2-26) 

T ' (x, Z , t) = T 1 (Z,t)cos yx + T 2 (Z,t)sin yx , (2-27) 
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(2-28) 



<j>^(x,t) = S 1 (t)cos yx + S 2 (t)sin yx , 


(2-28) 


u'(x,Z,t) = U^(Z,t)cos yx + U 2 (Z,t)sin yx , 


(2-29) 


v'(x,Z,t) = V^(Z,t)cos yx + V 2 (Z,t)sin yx , 


(2-30) 


Z'(x,2,t) = W^(Z,t)cos yx + W 2 (Z,t)sin yx , 


(2-31) 


<J>'(x,Z,t) = H^(Z,t)cos yx + H 2 (Z,t)sin yx , 


(2-32) 


Q'(x,Z,t) = Q 1 (Z,t)cos yx + Q 2 (Z,t)sin yx , 


(2-33) 


MTS'(x,t) = HTS^(t)cos yx + MTS 2 (t) sin yx . 


(2-34) 



The relations (2-25) - (2-34) are substituted into equations 



(2-15 )-( 2-22 ) . The prognostic and diagnostic equations 


are 


separated into equations for the cosine and sine terms, 
resultant prognostic equations are 


The 


3A 

3t = - (D 1 - “^2 - 8V 1 ' 


(2-35) 


3A, 

3t - - fD 2 + U|JA 1 * SV 2 • 


(2-36) 


3D. _ 2 

3t fA l uyD 2 ^ U 1 M dZ W 2 + V H 1 ' 


(2-37) 
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(2-38) 



3D 2 

3t" 



= fA 2 + uuD ± - 6U 2 + y g| W x + m 2 H 



9T] 

3t 



f du ^ 

-U ^ T 2 + R dZ V 1 - W 1 + Q 1 ' 



3T, 

W 



f r? 11 C} m 

" UtJT l + RdZ V 2 - 3Z W 2 + Q 2 ' 



9S ] 

at 



= -uys 2 + fuV 1 - RTW X + MTS^ , 



as, 

at 



uys x + fuV 2 - rtw 2 + MTS 2 



The resultant diagnostic equations for u and v are 



U 



U 



V 



V 



1 



2 



1 



2 




(2-39) 

(2-40) 

(2-41) 

(2-42) 

(2-43) 

(2-44) 

(2-45) 

(2-46) 



Geopotential values above the surface are obtained by inte- 
grating the hydrostatic equation from the surface (Z = Z^) to 
height Z: 

Z 

H, = R / T, (Z,t)dZ + S, (2-47) 

Z 0 
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(2-48) 



Z 

H 9 = R / T_ (Z, t)dZ + S 

z o 



The vertical velocity is calculated by integrating the con- 
tinuity equation from the top of the atmosphere (Z = Z T ) down 
to height Z. The upper boundary condition, Z = 0 at Z = Z T , 
is used. This boundary condition is not exact, but some form 
of it is used in most numerical models. The diagnostic 
equations for the vertical velocity are 



W x = / D x (Z,t)dZ , 



W 2 = / D 2 (Z,t)dZ . 



(2-49) 



(2-50) 



Equations (2-35) - (2-50) are the prognostic and diagnostic 
equations that govern all six numerical models. Using the given 
basic state and the one-wave spectral perturbation quantities, 
the governing equations reduce to functions of Z and t. The 
models are effectively one-dimensional. 

To display the results of each model, the sine and cosine 
amplitudes of each variable are combined to determine the 
amplitude and phase of a single cosine wave in the x-direction. 

A typical variable has the form 



Y (x , Z , t ) = A(Z, t) cos (yx-6) , 



(2-51) 
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where the amplitude is A(Z,t) and the phase is <5(Z,t). The 
amplitude and phase are calculated at each level for all 
variables . 

C. TIME DIFFERENCING 

Two forward time steps are taken to start each model and 
then leapfrog time differencing is used. The leapfrog scheme 
is employed because of its ease to code. A Robert filter is 
used to reduce the amplitude of the computational mode gener- 
ated by the leapfrog time differencing. The filter is dis- 
cussed by Haltiner and Williams, 1980. For a prognostic 
variable F, calculate F the average value of F at time 

step (n-l)At, using equation (2-52). Using the unaveraged 

OF 

values at time step nAt, compute the tendency (-s-p) from its 

dr n 

predictive equation. The predicted value at time step 
(n+l)At is calculated using equation (2-53). 



n-1 



F ,+y(F-2F -, + F _) 

n-1 ' n n-1 n-2 



(2-52) 



where y is a weighting function. 



n+1 



= F 



n-1 



+ 2At(||) 



n 



(2-53) 



In all thesis experiments, y = 0.05 is used. The time step 
for each experiment is calculated in the model by requiring, 
for computational stability. 
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vAt 



1 

2 



9 



( 2 - 54 ) 



where v = pc, and c is the typical phase speed of an external 
gravity wave. 

D. VERTICAL GRIDS 

Each of the models uses one of three vertical grids. The 
three ways of distributing the variables over discrete levels 
are depicted in Figure 1. The staggered levels are represented 
by the dashed lines in Figure 1. Notice that the heights at 
which the variables are defined change between the three 
grids. The notation used in this paper to denote the staggered 
and unstaggered levels is consistent with the conventions used 
in the coded models. The height of the unstaggered levels 
is denoted as Z'. The height of the staggered levels is 
denoted as Z. In the models, both Z| and Z^ are defined to 
be the surface of the earth. It is assumed that the staggered 
level Z^ is exactly in the middle of the layer between Z^_^ 
and Z|. This distinction is important because the models can 
have layers with unequal depth. Thus, the height of the 
staggered levels is defined relative to the height of the 
unstaggered levels. 

A finite difference model is written for each of the grid 
structures. The models are denoted as DFM-A , FDM-B and 
FDM-C. Similarly, finite element models using the three grids 
are indicated by FEM-A, FEM-B and FEM-C. 
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E. FINITE DIFFERENCE MODELS 

The only differences in the equations between the three 
FDM models are the approximations of terms involving du/dZ 
and 3T/3Z in the prognostic equations and the approximations 
of the integral in the diagnostic geopotential equation. 

Centered difference approximations are used, except at the 
boundaries where one-sided differences are employed. The 
finite difference approximations used in the prognostic equa- 
tions are listed in Appendix A. 

F. FINITE ELEMENT MODELS 
1. FEM-C 

The unstaggered FEM model is the simplest of the three 
FEM models. Each of the dependent variables is expanded into 
a finite series in terms of the eigenfunctions <|>j (Z) . The 
eigenfunctions for this model are depicted in Figure 2. The 
eigenfunction expansion for a typical term is 

N+l 

A-, ( Z , t ) = l A? ( t ) <j) . ( Z ) . (2-55) 

j=l 3 

The finite element approximations for the vorticity, 
divergence and thermodynamic equations are derived by substi- 
tuting the eigenfunction expansion for each dependent variable 
into equations (2-35) - (2-40) . Each equation is multiplied 
by 4> ± (Z) and integrated with respect to Z from the bottom to 
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Figure 2. Eigenfunctions defined at unstaggered levels 
(solid lines) for grid C. Dashed lines 
represent staggered levels. 



the top of the atmosphere. Each term in the equations is the 
finite sum of separate integrals. Only the integrals of 
overlapping eigenfunctions are nonzero. The resultant equa- 
tions, listed in Appendix B, are matrix equations. For an 
N-layer model, the variables are calculated at N+l discrete 
levels and the variable can be represented as a column vector 
of length N+l. Thus, the forecast equations for vorticity, 
divergence and potential temperature become vector equations 
which contain four (N+l) x (N+l) matrices multiplying appro- 
priate vectors. 
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The mass matrix for this model is defined by 



M = 



i +1 

l 

j=i-l 



/ <J> • (Z) <J> . (Z) dZ , for i = 1,...,N+1 



(2-56) 



The matrix N is defined for terms multiplied by u. 



i+1 i+1 



N = 



l l u. / 4> - ( Z ) ( Z ) 0 . (Z) dZ , 

j=i-l k=i-l 3 Z Q J k i 

for i = 1 , . . . , N+l . 

The matrix P is defined for terms multiplied by du/dZ , 



(2-57) 



i+1 - i+1 



J T d<f>. 



P = 



l l u. / dZ 1 <f> k (z)M z)dz ' 

j=i-l k=i-l 3 z Q 1 

for i = 1 , . . . , N+l . 

The matrix R is defined for terms multiplied by 9T/3Z, 



(2-58) 



R 



i+1 i+1 

l l T. 

j=i-l k=i-l J 



Z T d 4> • 

/ ^ 4> k (z) 4. i (z)dz, 



for i = 1, . . . ,N+1 . 



(2-59) 



The method to evaluate the terms of matrices ( 2-56 )-( 2-59 ) 
and the resultant formulas for the matrix elements are listed 
in Appendix C. 
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The vorticity, divergence and thermodynamic equations, 
written in matrix (indicated by a double underline) and 
vector (indicated by a single underline) form are 



M*|^(A2) 

(Dl) 



M* (~f *D1 - 6-V1) - y-N*A2 , 

M* (-f *D2 - B*V2) + yN-Al , 

M-(f-Al-BUl) ~ U*N*D2 - U-P-W2 
2 

+ y -M-H1 , 

M*(f*A2-BU2) + yN-Dl + u -P-W2 
+ y 2 *M* H2 , 



M~(T1) 

M*|r (T2) 



-y -N-T2 + (|) *P*V1 - R- W1 + M-Ql , 
y * N • Tl + (^) -P-V2 - R-W2 + M-Q2 . 



(2-60) 

(2-61) 

(2-62) 

(2-63) 

(2-64) 

(2-65) 



Equations (2-60) - (2-65) are simplified by multiplying 
each equation by M ^ and applying the Robert filter. The 
matrices M ^ * N, M ^*P, M ^*R are constants. They are con- 
structed in the initialization subroutine and stored for use 
in the forecast subroutine. The matrices are multiplied by 
the appropriate vectors with values for time level nAt. The 
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resultant forecast equations are vector equations and the fore- 
cast value for the i-th vertical level is the sum of values 
in the i-th location of each vector in the equation. The 
prognostic equations for the vorticity, divergence and poten- 
tial temperature vectors are 



-1 



Al n+1 - = AlBAR n 1 + 2 At (-f • D1 - B*V1 - y*M ^N*A2) 



( 2 - 66 ) 



-1 



A2 n+1 = A2BAR n _ 1 + 2 At ( -f • D2 - B • V2 + y • M ♦N * Al ) 



(2-67) 



D1 , •» 
— n+1 



-1 



= DlBAR n _ 1 + 2At(f -Al - B'Ul " y’M -N-D2 



- y • M 1 • P • W2 + y 2 Hl) , 
- n 



( 2 - 68 )' 






D2BAR . + 2 At ( f • A2 - B * U2 + y • M~ 1 • N • Dl 
n-1 — — — — 



+ y*M 1 -P-W1 + y 2 H2) 

~~ — — n 



(2-69) 



—n+1 



T2 

— n+1 



TlBAR . + 2At (-y •M~ 1 *N-T2 
n-1 — — 

+ (|) • m" 1 ♦ P • VI - M _1 • R * W1 + Ql) , 



T2BAR . + 2 At ( y * M 1 *N*T1 
n-1 - — 

+ (1)*M -1 *P*V2 - M _1 -R-W2 + Q2) , 

R — — - n 



(2-70) 



(2-71) 



where the subscripts n+1, n, n-1, refer to the values of the 
vectors at time step (n+1) At, nAt, and (n-1) At, respectively. 
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The surface geopotential and the diagnostic variables 
are calculated using the same equations that are used in 
the unstaggered model FDM-C. 

2. FEM-A 

The FEM-A model defines vertical velocity at the un- 
staggered levels in terms of the eigenfunctions \p^(Z). The 
other variables are defined at the staggered levels in terms 
of the eigenfunctions 4>j (Z) . The eigenfunctions for this 
model are depicted in Figure 3. 




Figure 3. Eigenfunctions for grids A and B. Egienf unc- 
tions <MZ) are defined for the unstaggered 
levels (solid lines at height Z') and $(Z) 
are defined for the staggered levels 
(dashed lines at height Z) . 
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The finite element approximations for the vorticity, 
divergence and thermodynamic equations are derived using the 
method described for the FEM-C model. The resultant Galerkin 
formulation of the equations, listed in Appendix D, are matrix 
equations. The mass matrix, M, and the matrix, N, for terms 
of the form u multiplied by variables defined at the nodal 
points of <f>j (Z) have the same formulas as M and N in the FEM-C 
model, equations (2-56) and (2-57) . However, the eigenfunc- 
tions 4>j (Z) are not defined at the same levels in both models. 

dii 

The matrix Q is defined for terms of the form 



Q 



i+1 i+1 

l l Uj 

j=i-l k=i-l 3 



Z T d<*> . 

/ dZ 1 Mi dZ 



(2-72) 



The matrix P is defined for terms of the form 3 -^*W, 
= dZ 



P = 



i+1 • i+1 _ T d4> . 

I l u i / dT 1 Mi dz 

j=i-l k=i-2 3 Z n K 1 



(2-73) 



9T 



The matrix R is defined for terms of the form 



i+1 i+1 



J T dcf> . 



R = 



l - - i , T j „ / dz Vi dz 



(2-74) 



j=i-l k=i-2 J Z. 



The staggered eigenfunctions present three general 
problems for evaluating the elements of the five matrices. 
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First, for an N-layer model, portions of eigenfunctions 
4> 1 (Z) and (^^(Z) are defined in the model atmosphere but 
the physical meaning of contributions from those terms is 
unclear. ' The contributions are included in the first two 
rows and the last row of each matrix. Second, only portions 
of eigenfunctions (^(Z) and 4>^ + ^(Z) are defined in the model 
atmosphere. To describe the incomplete sides of both eigen- 
functions an assumption must be made about the value of 4> ^ 
at the surface and <t> N+ ^ at the top of the atmosphere. Last, 
the equations derived for the general elements of each matrix 
are more complex than the formulas for the model FEM-C. The 
complexity of these equations makes it much more difficult 
to evaluate the elements in the first two and last two rows 
of each matrix. 

Assumptions are made and procedures are developed in 
an attempt to resolve these problems. In this model the mean 
state variables, u and T, are defined only at the N staggered 
levels. However, u and T values defined at the nodal points 
of 0^(Z) and (2) are important in the Galerkin formulation 

of the du/dZ and 9T/9Z terms. In my thesis experiments, the 
values of u and T are defined at the surface and top of the 
atmosphere; they are not defined at the nodal points of <f>^(Z) 
and 4 ) n+ 2 ^ z ^ becuase these points are outside the model atmos- 
phere. For constant shear with height, u and T are defined 
at the boundaries such that the shear in the two half layers 
at the boundaries is the same as the shear in the other layers. 
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The contributions from the perturbation quantities defined 
at the nodal points of <j> 1 (Z) and < J , N+ 2 ^ Z ^ are not included in 
these experiments. To be consistent with the other N-layer 
models, all perturbation quantities are defined as column 
vectors of length N+l even though all quantities, except ver- 
tical velocity, are defined at only N levels. The values of 
each perturbation at its lowest staggered level are stored 
as the second element in its column vector; zero is stored 
as the first element, which corresponds to the contribution 
from <j>^(Z). To evaluate the staggered eigenfunctions defined 
in the layers between the surface and an< 3 z n+ 1 an< ^ 

top of the atmosphere, it is assumed that the value at the 
boundaries of those eigenfunctions is one-half. Thus, three- 
fourths of the eigenfunctions ^ ( 2 ) and <i> N+ 2^ z ^ are defined 
in the model atmosphere. 

The equations for the general elements of the five 
matrices are evaluated by substituting into equations (2-56)- 
(2-57) and (2-72) - (2-74 ) the formulas for ({^^(Z), <J>^(Z), 

4>i_i (Z ) , * (Z), ^(Z), * (Z), and <Jk_ 2 (Z) terms of the 

local coordinate £ = Z -Z^. The equations for these eigen- 
functions defined for the levels 1, 2, i, and N+l are listed 
in Appendix E. The equations for the elements of the first 
two rows and the last row of the matrices P, Q, and R are 
derived using the formulas for levels 1, 2, and N+l. However, 
the corresponding rows of matrices M and N were not calcu- 
lated using these formulas because the model results for the 
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barotropic experiment were not constant with height when the 
elements in those three rows of both M and N were calculated 
using these formulas. The reason for this is not clear. The 
equations for the general elements of M and N are simple 
enough to allow the terms in the affected three rows of each 
matrix to be deduced by reason. The equations for the elements 
of all five matrices are listed in Appendix F. 

The prognostic matrix equations for vorticity, diver- 
gence and temperature, listed in Appendix D, are simplified 
in a manner similar to the method described for the- FEM-C 
model. The final form of the forecast equations for the 
vorticity and divergence vectors are the same as for model 
FEM-C, equations (2-66) - (2-69) . The thermodynamic vector 
equations are 

TlBAR n _i + 2At ♦ [Q1 - p- M~ 1 -N - T2 

+ (J) •M _ 1 -Q-V1 - M - 1 *R-W1] , (2-75) 

K ■ ■ - — ~ — — n 

T2BAR n + 2At • [Q2 + U' M" 1 -N - T1 

+ (- 5 -) M _1 * Q • V2 - M _1 • R* W2] . (2-76) 

R ■ = = — — n 

The six forecast equations are solved at levels two to N+l . 

The surface geopotential and the diagnostic variables are 
calculated using the corresponding equations in model FEM-A. 



T1 .. 
— n+l 



T2 i n 
— n+l 
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3. FEM-B 



The FEM-B model defines vertical velocity, potential 



temperature, mean state potential temperature and diabatic 
heating at the unstaggered levels in terms of the eigenfunc- 



tions are the same as defined for the FEM-A model, shown in 
Figure 3. 



divergence and thermodynamic equations are. derived by substi- 
tuting the eigenfunction expansion for each dependent variable 
into equations (2-35) - (2-40) . The vorticity and divergence 
equations are multiplied by <J> ^ ( Z ) and integrated with respect 
to Z from the bottom to the top of the atmosphere. The resul- 
tant Galerkin formulation of the vorticity and divergence 
equations are the same as those derived for model FEM-A. The 
matrices in those equations, M, N, and P, are the same as 
defined for FEM-A, equations (2-56) - (2-57) and (2-73). The 
thermodynamic equations are multiplied by i(^(Z) because poten- 
tial temperature is defined at the unstaggered levels. As 
before, the equations are integrated through the depth of the 
atmosphere. The resultant equations are listed in Appendix 
G. Four additional matrices are defined for the two thermo- 
dynamic equations. The mass matrix S is 



tions (Z) . The other variables are defined at the staggered 
levels in terms of the eigenfunctions j(Z). The eigenfunc- 



The finite element approximations for the vorticity 



i+1 



S 




/ j (Z) (Z) dZ . 



(2-77) 
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The matrix Q is defined for terms of the form ^r’V, 



i+2 i+2 



Q = 



I I u. / <j>. (Z) <*> k (Z) (Z) 

i=i-l k=i-l J Z. J 



dZ 



( 2 - 78 ) 



B T 

The matrix R is defined for terms of the form -^=-‘W, 

— 0 L 



i+1 i+1 



J T di{>_ 



R = I It, / ^^(zi^izjdz 

j=i-l k=i-l 3 Z Q dZ k 1 



( 2 - 79 ) 



The matrix T is defined for terms of the form u*T, 



i+2 i+1 



T = 



I I u. J 4) • (Z)^ k (Z)^.(Z) 

i=i-l k=i-l J Z n 3 



dZ 



( 2 - 80 ) 



As discussed in the FEM-A model description, the 
staggered finite elements present problems for evaluating the 
elements of the matrices. In this model u is defined at the 
surface, the top of the atmosphere, and at the N staggered 
levels. The mean state temperature, T, is defined at the 
unstaggered levels so special definitions for it are not needed. 
Also, the contributions from the perturbation quantities 
defined at the nodal points of 4>^(Z) and are not ^ n- 

cluded in the model FEM-B model experiments. The staggered 
eigenfunctions, cJ)^(Z), are evaluated at the boundaries using 
the assumptions discussed in the previous section. 
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The elements of matrices S, Q, R, and T are evaluated 
by substituting formulas for <{> i+2 • 4> i+1 (Z), <J>j.(Z), 

<|k +1 (Z), i/j ^ ( Z ) and ^ i _ 1 (Z) / defined in terms of the local 
coordinate n = Z -Z^, into equations (2-77 )-( 2-80 ) . Formulas 
for these eigenfunctions are listed in Appendix H. The equa- 
tions for the matrix elements are listed in Appendix I. 

The forecast matrix equations for vorticity, divergence 
and temperature are simplified in a manner similar to the 
method described for model FEM-C. The final form of the vor- 
ticity and divergence vector equations are the same as for 
model FEM-C, equations (2-66) - (2-69) . The thermodynamic vec- 
tor equations are 

TlBAR n + 2 At [Q1 - y S~ 1 -T -T2 

+ (i)-S _1 -Q-Vl - S _1 • R*W1] , (2-81) 

k - ~ == — n 

T2BAR n + 2 At [Q2 + y S~ 1 T - Tl 

+ (J)-S -1 -Q-V2 - S _1 * R* W2 ] . (2-82) 

k . ,, ■ — 1 " 1 — n 

The vorticity and divergence equations are solved at levels 
two to N+l . The thermodynamic equations are solved at all 
unstaggered levels. The surface geopotential and the diagnos- 
tic variables are calculated using the corresponding equations 
in model FDM-B. 



T1 .li 

— n+l 



T2 . . 
— n+l 



38 



III. EXPERIMENTS AND RESULTS 



Three experiments are performed with each model; an 
initial perturbation in the meridional flow, flow over moun- 
tain topography and flow with a diabatic heat source. The 
first two experiments are performed with each model defined 
with six and then sixty layers. The heating experiment is 
repeated with six, twelve and sixty-layer models. The analytic 
solution of each experiment has not been derived. For each 
experiment, the 60-layer model results are intercompared to 
determine if the models are converging to the same solution. 

The standard of comparison for each six-layer model is its 
corresponding 60-layer solution. Temperature and divergence 
profiles are examined in each experiment. 

Several parameters are defined identically in each experi- 
ment. The vertical coordinate, Z, is defined between zero 
and one (1000-368 mb) and the vertical levels are equally 
spaced. The x-wavelength is 4,000 kilometers. The time step 
is 17.7 minutes. The Coriolis parameter is defined at 45 
degrees latitude. There is no vertical shear in the u field 
and u = 10.0 meters/second (m/s). The mean state potential 
temperature increases with height from its surface value of 
310.0 Kelvin (K) . 

A. ROSSBY WAVE EXPERIMENT 

Rossby waves are generated in each model using an initial 
perturbation, v' = 5.0 m/ s , in the cosine term of the 
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meridional flow. All other perturbations are initially zero. 
There is no diabatic heat source and no mountain topography. 

The latitudinal variation of the Coriolis parameter, 6, is 
defined at 45 degrees latitude. The forecast length of this 
experiment is 96 hours. 

1 . Sixty-Layer Models 

"The 60-layer FDM-A, FDM-B, FDM-C, and FEM-C models 
converge to the same temperature and divergence solutions 
(Figures 4-7) . It should be noted that the phase of each 
variable is defined between zero and 360 degrees. There is 
a discontinuity in the phase profile if the phase passes 
through zero degrees. The height at which the temperature 
phase discontinuity in model FDM-A occurs differs from the 
other three models because temperature is defined at the 
staggered levels in FDM-A. The four models represent the 
same physical solution, which is called the consensus solution. 

The FEM-B temperature amplitude is slightly larger 
than the consensus amplitude (5%) and an amplitude oscillation 
is present in the lowest three layers of the atmosphere (Figure 
8) . The jagged profile may be caused by the terms in the 
matrices which represent contributions from the eigenfunctions 
near the lower boundary of the model. There is a discontinuity 
between the temperature phases in the lowest two layers of the 
FEM-B profile (Figure 9). Apart from this feature, the 
general shape of the phase profile is similar to the consensus, 
although the phase is five degrees less and passes through 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (60 LAYERS) 




LEGEND 
□ - FDM-A 
o = FDM-B 

* = FDM-C 

♦ = FEM-C 



Figure 4. Sixty-layer Rossby-wave experiment at 96 

hours. Temperature amplitude profiles are 
compared for models FDM-A, FDM-B, FDM-C, 
and FEM-C. 
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HEIGHT. Z = - LN(P/PO) 



TEMPERATURE PHASE FOR V = 5.0 CASE (60 LAYERS) 




TEMPERATURE PHASE (0-360 DEGREES) 



LEGEND 
c = FDM-A 
c = FDM-B 
£ = FDM-C 
* = FEM-C 



Figure 5. As in Fig. 4 but for temperature phase 



HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR V = 5.0 CASE (60 LAYERS) 




LEGEND 
□ _ FDM-A 
c = FDM-B 
a = FDM-C 
+ = FEM-C 



Figure 6. As in Fig. 4 but for divergence amplitude. 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE PHASE FOR V = 5.0 CASE (60 LAYERS) 




LEGEND 
□ = FDM-A 
c = FDM-B 
a = FDM-C 
* = FEM-C 



Figure 7. As in Fig. 4 but for divergence phase. 
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HEIGHT. Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (60 LAYERS) 




LEGEND 
= = FEM-B 
o = FDM-C 



Figure 8. Sixty-layer Rossby wave experiment at 96 

hours. FEM-B temperature amplitude profile 
is compared with the temperature amplitude 
profile of FDM-C, which represents the 
consensus profile. 
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HEIGHT. Z = - LN(P/PO) 



TEMPERATURE PHASE FOR V = 5.0 CASE (60 LAYERS) 




° = FDM-C 



Figure 9. As in Fig. 8 but for temperature phase. 
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zero at a lower level than the consensus. The divergence 
amplitude of model FEM-B is within five percent of the con- 
sensus (Figure 10) but its phase profile (Figure 11) is quite 
different from the consensus. 

The model FEM-A temperature and divergence ampli- 
tudes vary significantly from the consensus. The temperature 
amplitude profile (Figure 12) has a large amplitude oscilla- 
tion in the lowest three levels above the surface, possibly 
caused by the choice of eigenfunction representations used 
to evaluate elements of the first two rows of the matrices. 

The FEM-A temperature amplitude is larger than the consensus 
throughout the atmosphere. The temperature phase passes 
through zero at a lower level and is five degrees less than 
the consensus (Figure 13). The maximum FEM-A divergence 
amplitude occurs near the middle of the atmosphere and its 
magnitude is ten percent less than the maximum of the consensus 
(Figure 14). The consensus divergence phase is nearly con- 
stant with height, but the phase profile for model FEM-A 
steadily decreases with height (Figure 15) . 

2 . Six-Layer Models 

The comparison of six and sixty-layer profiles for 
variables defined at staggered levels may be initially mis- 
leading. The first staggered level in a six-layer model occurs 
at Z = 0.0833. The lowest staggered level in a 60-layer 
model is defined at Z = 0.0083, which may be mistaken for the 
surface in the graphs. When the values of a six-layer model 
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HEIGHT, 2 = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR V = 5.0 CASE (60 LAYERS) 




10.0 12.5 15.0 17.5 20. C 22.5 25.0 27.5 

DIVERGENCE AMPLITUDE (1/SEC) *10 9 



LEGEND 
r. = FEM-B 
o = FDM-C 



Figure 10. As in Fig. 8 but for divergence amplitude. 



HEIGHT. Z = - LN(P/PO) 



DIVERGENCE PHASE FOR V = 5.0 CASE (60 LAYERS) 




LEGEND 
□ = FEM-B 
c = FDM-C 



Figure 11. As in Fig. 8 but for divergence phase. 
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HEIGHT. Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (60 LAYERS) 




LEGEND 
□ = FEM-A 
c = FDM-C 



Figure 12. As in Fig. 8 but for model FEM-A compared 
with the consensus, FDM-C. 
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HEIGHT. Z = - LN(P/PO) 



TEMPERATURE PHASE FOR V = 5.0 CASE (60 LAYERS) 




Figure 13. As in Fig. 12 but for temperature phase. 
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HEIGHT. Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR V = 5.0 CASE (60 LAYERS) 




LEGEND 
a = FDM-C 
o = FEM-A 



Figure 14. As in Fig. 12 but for divergence amplitude. 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE PHASE FOR V = 5.0 CASE (60 LAYERS) 




LEGEND 
n = FEM-A 
c = FDM-C 



Figure 15. As in Fig. 12 but for divergence phase 
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coincide with a 60-layer profile, the models are considered 
to represent the same physical solution, even though the 
six-layer model has a smaller vertical domain for staggered 
variables. 

The temperature amplitude profiles of the six-layer 
FDM-A, FDM-B, FDM-C and FEM-C models are identical with their 
corresponding 60-layer results (Figures 16-19) . The previously 
discussed problems in the lowest layers of models FEM-A and 
FEM-B are quite evident in their six-layer profiles. The 
temperature amplitude of FEM-B is virtually identical with 
its 60-layer profile above the first two layers (Figure 20) . 

The six-layer FEM-A profile (Figure 21) does not agree with 
its 60-layer profile below Z = 0.50 (600 mb). The large 
amplitude oscillation in this profile destroys the integrity 
of the solution in a significant portion of the atmosphere. 

The six-layer divergence amplitude profiles for the 
grid C models are identical with each other (Figures 22-23) , 
and the divergence profiles for the grid B models are nearly 
the same (Figures 24-25). All four models approximate their 
corresponding 60-layer solutions with similar accuracy. The 
six-layer FDM-A model does not approximate the curvature of 
the consensus 60-layer divergence profile as well as the 
grid B and C models (Figure 26). The six-layer FEM-A model 
has difficulty approximating the location and magnitude of 
the mid-atmosphere divergence maximum of its 60-layer solution 
(Figure 27) . 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (FDM-A) 




LEGEND 

- = FDM-A 60 LAYERS 
c = FDM-A 6 LAYERS 



Figure 16. Six-layer Rossby wave experiment at 96 

hours. Temperature amplitude profiles are 
compared for the six-layer and 60-layer 
FDM-A models. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (FDM-B) 




LEGEND 

c = FDM-B 60 LAYERS 
o = FDM-B 6 LAYERS 



Figure 17. As in Fig. 16 but for model FDM-B. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (FDM-C) 




LEGEND 

= = FDM-C 6 LAYERS 
o = FDM-C 60 LAYERS 



Figure 18. 



As in Fig. 16 but for model FDM-C. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (FEM-C) 




LEGEND 

= = FEM-C 60 LAYERS 
c = FEM-C 6 LAYERS 



Figure 19. As in Fig. 16 but for model FEM-C. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (FEM-B) 




LEGEND 

c = FEM-B 60 LAYERS 
c = FEM-B 12 LAYERS 



Figure 20. As in Fig. 16 but for model FEM-B. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR V = 5.0 CASE (FEM-A) 




LEGEND 

c — FEM— A 60 LAYERS 
c = FEM— A 6 LAYERS 



Figure 21. As in Fig. 16 but for model FEM-A. 



60 



HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR V = 5.0 CASE (FDM-C) 




LEGEND 

c = FDM-C 60 LAYERS 
c = FDM-C 6 LAYERS 



Figure 22. As in Fig. 16 but divergence amplitude of 
FDM-C. 
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Figure 23. 



As in Fig. 22 but for model FEM-C. 



62 



HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE EOR V = 5.0 CASE (FDM-B) 




LEGEND 

□ = FDM-B 60 LAYERS 
c = FDM-B 6 LAYERS 



Figure 24. As in Fig. 22 but for model FDM-B. 
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DIVERGENCE AMPLITUDE FOR V = 5.0 CASE (FEM-B) 
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Figure 25. As in Fig. 22 but for model FEM-B. 
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Figure 26. As in Fig. 22 but for model FDM-A. 



65 



HEIGHT. Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR V = 5.0 CASE (FEM-A) 




LEGEND 

c = FEM-A 60 LAYERS 
o = FEM-A 6 LAYERS 



Figure 27. As in Fig. 22 but for model FEM-A. 
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B. MOUNTAIN TOPOGRAPHY EXPERIMENT 



The forced vertical velocity term, MTS, in the surface 
geopotential forecast equation (2-3) is non-zero in this 
experiment. It represents the contribution to surface geo- 
potential from air flowing over mountain topography which 
varies sinusoidally in the x-direction and has no variation 
in the y-direction. The mountain ridge-to-valley height 
difference is 1,500 meters. To reduce the trauma for the 
model, the mountains are gradually "built" to their full 
height over a period of 36 hours. Thus, the forced vertical 
velocity increases in the first 36 hours of the forecast 
period and is constant for the remainder of the 96-hour fore- 
cast period. The equations used to define the forced verti- 
cal velocity are included in Appendix J. Beta and all initial 
perturbations are zero in this experiment. 

1 . Sixty-Layer Models 

The FDM-A, FDM-B, FDM-C, and FEM-C models converge 
to the same physical solution. The amplitude and phase pro- 
files of temperature and divergence are depicted in Figures 
28-31. 

The FEM-B model again has a jagged temperature ampli- 
tude profile in the lowest two layers (Figure 32). The 
temperature amplitudes in the remainder of the atmosphere are 
within 5% of the consensus. The FEM-B temperature phase has 
oscillations in the bottom three layers and the top two layers 
of the profile (Figure 33) , but the rest of the profile is 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR MOUNTAIN CASE (60 LAYERS) 




TEMPERATURE AMPLITUDE (KELVIN) 



LEGEND 
□ = FDM-A 
o = FDM-B 
a = FDM-C 
* = FEM-C 



Figure 28. Sixty-layer mountain topography experiment 

at 96 hours. Temperature amplitude profiles 
are compared for models FDM-A, FDM-B, FDM-C, 
and FEM-C. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE PHASE FOR MOUNTAIN CASE (60 LAYERS) 




LEGEND 
c = FDM-A 
o = FDM-9 
a = FDM-C 
* = FEM-C 



Figure 29. As in Fig. 28 but for temperature phase. 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR MOUNTAIN CASE (60 LAYERS) 




LEGEND 
o = FDM-A 
c = FDM-B 
a = FDM-C 
* = FEM-C 



Figure 30. As in Fig. 28 but for divergence amplitude. 
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HEIGHT. Z = - LN(P/PO) 



DIVERGENCE PHASE FOR MOUNTAIN CASE (60 LAYERS) 
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a = FDM-C 
* = FEM-C 



Figure 31. As in Fig. 28 but for divergence phase. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR MOUNTAIN CASE (60 LAYERS) 




LEGEND 
c = FEM — B 
o = FDM-C 



Figure 32. Sixty-layer mountain topography experiment 
at 96 hours. FEM-B temperature amplitude 
profile is compared with the temperature 
amplitude profile of FDM-C, which represents 
the consensus profile. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE PHASE FOR MOUNTAIN CASE (60 LAYERS) 




TEMPERATURE PHASE (0-360 DEGREES) 



LEGEND 
a = FEM-B 
o = FDM-C 



Figure 33. As in Fig. 32 but for temperature pha 



within one degree of the consensus phase. The FEM-B 
divergence amplitude and phase are identical with the 
consensus values. 

An oscillation also exists in the lowest three layers 
of the temperature amplitude profile for model FEM-A (Figure 
34) . The amplitude is 30% higher than the consensus near the 
surf ace, and this difference decreases with height. The tem- 
perature phases are nearly identical. The FEM-A divergence 
amplitude profile is similar to the consensus, but the amplitude 
does not decrease as fast with height as the consensus 
(Figure 35) . The divergence phase profile for this model is 
the same as the consensus. 

2 . Six-Layer Models 

The temperature amplitude profiles of the six-layer 
models FDM-A, FDM-B, FDM-C and FEM-C are identical with each 
other and also with the consensus solution (Figures 36-39) . 

The six-layer staggered finite element models are nearly 
identical with their 60-layer solutions above the first two 
layers (Figures 40-41) . The previously identified inadequate 
representations of the temperature amplitude in the lowest 
layers exist in this experiment. All six models approximate 
their 60-layer divergence amplitude profiles quite well 
(Figures 42-46) , however the FEM-A model (Figure 47) does 
not approximate the curvature in the lower portion of its 
60-layer profile as well as the other models. 
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HEIGHT, 2 = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR MOUNTAIN CASE (60 LAYERS) 




LEGEND 
□ = FEM-A 
e = FDM-C 



Figure 34. As in Fig. 32 but for model FEM-A compared 
with FDM-C. 
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HEIGHT. Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR MOUNTAIN CASE (60 LAYERS) 




LEGEND 
c = FEM-A 
o = FDM-C 



Figure 35. As in Fig. 34 but for divergence amplitude. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR MOUNTAIN CASE (FDM-A) 




LEGEND 

a = FDM-A 60 LAYERS 
o = FDM-A 6 LAYERS 



Figure 36. Six-layer mountain topography experiment at 
96 hours. Temperature amplitude profiles 
are compared for the six-layer and 60-layer 
FDM-A models. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR MOUNTAIN CASE (FDM-B) 




TEMPERATURE AMPLITUDE (KELVIN) 



LEGEND 

= = FDM-3 60 LAYERS 
o = FDM-3 6 LAYERS 



Figure 37. As in Fig. 36 but for model FDM-B. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR MOUNTAIN CASE (FDM-C) 




TEMPERATURE AMPLITUDE (KELVIN) 



LEGEND 

= FDM-C 6 LAYERS 
c = FDM-C 60 LAYERS 



Figure 38. As in Fig. 36 but for model FDM-C 
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Figure 39. As in Fig. 36 but for model FEM-C 



HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR MOUNTAIN CASE (FEM-A) 




LEGEND 

c = FEM-A 60 LAYERS 
o = FEM-A 6 LA'tERS 



Figure 40. As in Fig. 36 but for model FEM-A. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR MOUNTAIN CASE (FEM-B) 
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c = EEM-B 60 LAYERS 
o = FEM-B 6 LAYERS 



Figure 41. 



As in Fig. 36 but for model FEM-B 



HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR MOUNTAIN CASE (FDM-A) 




LEGEND 

n = FDM-A 60 LA'i ERS 
c = FDM-A 6 LAYERS 



Figure 42. As in Fig. 36 but for divergence amplitude. 
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DIVERGENCE AMPLITUDE TOR MOUNTAIN CASE (FDM-B) 




LEGEND 

c = FDM-B 60 LAYERS 
c = FDM-B 6 LAYERS 



Figure 43. As in Fig. 36 but divergence amplitude 
for model FDM-B. 
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Figure 44. 



As in Fig. 36 but divergence amplitude 
for model FDM-C. 
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DIVERGENCE AMPLITUDE FOR MOUNTAIN CASE (FEM-B) 
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Figure 45. As in Fig. 36 but for divergence amplitude 
for model FEM-B. 
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DIVERGENCE AMPLITUDE FOR MOUNTAIN CASE (FEM-C) 
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o = FEM-C 6 LAYERS 



Figure 46. As in Fig. 36 but for divergence amplitude 
for model FEM-C. 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR MOUNTAIN CASE (FEM-A) 




LEGEND 

c = FEM-A 60 LAYERS 
c = FEM-A 6 LAYERS 



Figure 47. 



As in Fig. 36 but for divergence amplitude 
for model FEM-A. 
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C. DIABATIC HEATING EXPERIMENT 



A diabatic heat source is defined in the layer between 
Z = 0.40 and Z = 0.60 (670-549 mb). The rate of heating is 
constant. in time and varies in x and Z, 



Q (x, Z , t) 



HEATING* cos 2 



(z -V 

It 7 7 — J * TT J -cos (yx) , 

(Z U Z L) 



(3-1) 



where HEATING is 5.0 K/day, Z^ is the midpoint of the heated 
layer, and Z L and Z y are the lower and upper boundaries, 
respectively, of the heated layer. The diabatic heating 
vectors, Q1 and Q2_, are defined in the initialization subrou- 
tine and stored for use in the forecast subroutine. Beta and 
all initial perturbations are zero. The forecast length is 
12 hours. 

1 . Sixty-Layer Models 

For the diabatic heating function defined in equation 
(3-1), the maximum heating occurs at Z = 0.50, the midpoint 
of the heated layer. The models defined using grids B and C 
define temperature and the heating functions at this point. 

The grid A models do not have temperature and diabatic heating 
defined at this point so the maximum rate of heating in these 
models is slightly less than in the other models, and the 
maximum heating occurs throughout one layer rather than 
occurring at one point. The heating rate at each level is 
listed in Appendix K for the six, twelve and sixty-layer 
staggered and unstaggered models. 
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The 60-layer profiles for the six models are quite 
similar, but the differences occur because the models are 
responding to different forcing. The temperature amplitude 
profiles -for the B and C grids come to a sharp point at 
Z = 0.50 and the grid A models have a square-nosed profile 
around this point (Figures 48-50) . The height in the atmos- 
phere ai; which the model experiences the onset of forcing is 
slightly different between the grids. This explains the small- 
scale oscillation in the temperature amplitude profiles near 
the boundaries of the heated layer, Z = 0.40 and Z = 0.60. 

The previously identified temperature amplitude oscillations 
in the lowest layers of models FEM-A and FEM-B are not evident 
in this experiment. The temperature phase profiles are nearly 
identical, except that the grid A models do not have as sharp 
of a spike at the boundaries of the heated layer as the other 
grids (Figures 51-53) . In summary, the 60-layer temperature 
profiles of all six models represent the same physical response 
to the diabatic heating. 

Divergence is defined at the midpoint of the heated 
layer only in the grid C models. Consequently, the shape of 
the divergence amplitude profile for the grid C models is 
different than the other four models. The grid C profiles 
have a sharp point at Z = 0.50 (Figure 54) . However, the 
minimum divergence is not symmetric around Z = 0.50 in the 
grid A and B models (Figures 55-57). The divergence ampli- 
tude is identical outside the heated layer for all models 
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Figure 48. Sixty-layer diabatic heating experiment at 
12 hours. Temperature amplitude profiles 
are compared for models FDM-A, FDM-B, FDM-C, 
and FEM-C. 
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Figure 49. Sixty-layer diabatic heating experiment at 

12 hours. FEM-A temperature amplitude profile 
is compared with the temperature amplitude 
profile of FDM-C, which represents the 
consensus profile. 
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pj_gu £“0 50. As in Fig. 49 but fox’ model FEM B compmred 
with the consensus, FDM-C. 
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TEMPERATURE PHASE FOR HEATING CASE (60 LAYERS) 
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* = FEM-C 



Figure 51. As in Fig. 48 but for temperature phase. 
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Figure 52. As in Fig. 49 but for temperature phas 



TEMPERATURE PHASE FOR HEATING CASE (60 LAYERS) 




LEGEND 
c = FEM — B 
o = FDM-C 



Figure 53. As in Fig. 50 but for temperature phase. 
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DIVERGENCE AMPLITUDE FOR HEATING CASE (60 LAYERS) 
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Figure 54. Sixty-layer diabatic heating experiment at 
12 hours. Temperature amplitude models are 
compared for the two grid C models. 
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Figure 55. Sixty-layer diabatic heating experiment at 
12 hours. Temperature amplitude profiles 
are compared for models FDM-A and FDM-B. 
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DIVERGENCE AMPLITUDE FOR HEATING CASE (60 LAYERS) 
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Figure 56. As in Fig. 55 but models FEM-B and FDM-C 
are compared. 
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Figure 57. As in Fig. 55 but models FEM-A and FDM-C 
are compared. 
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except FEM-A (Figure 57) . The divergence phase profiles are 
virtually identical for models FDM-A, FDM-B , FDM-C and FEM-C 
(Figure 58) . The divergence phase passes through zero at a 
different level for model FEM-B than the consensus, but this 
difference is not physically significant (Figure 59). The 
FEM-A profile differs slightly from the consensus outside the 
heated ^ayer (Figure 60). In conclusion, the six models are 
converging to the same 60-layer divergence solution inside 
the heated layer, but model FEM-A does not agree with the 
solution of the other five models outside this region. 

2 . Six and Twelve-Layer Models 

The difference between grids is more evident in this 
experiment than in the other experiments. Each model is run 
with both six and twelve layers and the results are compared 
with the 60-layer consensus of the six models. The tempera- 
ture amplitude consensus profile has a small-scale decrease 
at each boundary of the heated layer that cannot be reproduced 
using either six- or twelve-layer resolution. 

The six-layer grid A model temperature and divergence 
fields barely respond to the diabatic heating. The grid A 
models have identical temperature responses for both six and 
twelve layers. The six-layer perturbation, constant with 
height in the heated layer, is an order of magnitude smaller 
than maximum consensus perturbation (Figures 61-62). The 
grid A models have a stronger response to the heating with 
twelve layers than with six layers, but the maximum amplitude 
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Figure 58. As in Fig. 48 but for divergence phase. 
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Figure 59. As in Fig. 56 but for divergence phase. 
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DIVERGENCE PHASE FOR HEATING CASE (60 LAYERS) 
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Figure 60. As in Fig. 57 but for divergence phase. 
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Figure 61. Six-layer diabatic heating experiment at 12 
hours. Temperature amplitude profiles are 
compared for the six-layer and 60-layer 
FDM-A models. 
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Figure 62. As in Fig. 61 but for model FEM-A. 
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is only 60% of the 60-layer consensus response (Figures 
63-64). For both the six- and twelve-layer profiles, the 
temperature amplitude at the two staggered nodal points 
defined closest to the heated layer lie directly on the 60- 
layer profile. This may indicate that for much less than 
60-layer resolution the amplitude of the spike will not equal 
the maximum amplitude in the 60-layer experiment. The 
divergence fields of the six-layer grid A models is also an 
order of magnitude lower than their corresponding 60-layer 
solutions (Figures 65-66) . The characteristic sharp divergence 
decrease in the heated layer is roughly approximated by the 
12-layer models (Figures 67-68) . The minimum divergence is 
much larger than in the consensus profile and the base of the 
spike decrease is much broader than in the 60-layer solutions. 
Outside the heated layer, the shape of the 12-layer profile 
is similar to the consensus, but the amplitude is lower. 

The six-layer temperature amplitude response of the 
B and C grid models are virtually identical (Figures 69-72) . 

The base of the profile spike is broader in the six-layer 
models than in the consensus, but the four models closely 
approximate the magnitude and width of the tip of the spike. 

The 12-layer temperature profiles for grids B and C are also 
identical (Figures 73-76) . The profiles slightly underesti- 
mate the magnitude of the perturbation in the top half of the 
spikes, but the width at the base of the spikes are the same 
as in the 60-layer consensus. The small scale dip at the 
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TEMPERATURE AMPLITUDE (KELVIN) 



LEGEND 

c = FDM-A 60 LAYERS 
c = FDM-A 12 LAYERS 



Figure 63. Twelve-layer diabatic heating experiment at 
12 hours. Temperature amplitude profiles 
are compared for the 12-layer and 60-layer 
FDM-A models. 
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HEIGHT, 2 = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR HEATING CASE (FEM-A) 




LEGEND 

c = FEM-A 60 LAYERS 
o = FEM-A 12 LAYERS 



Figure 64. As in Fig. 63 but for model FEM-A 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR HEATING CASE (FDM-A) 




LEGEND 

= = FDM-A 60 LAYERS 
c = FDM-A 6 LAYERS 



Figure 65. As in Fig. 61 but for divergence amplitude. 



HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR HEATING CASE (FEM-A) 




0 T 

2.5 5.0 7.5 10.0 12.5 15.0 *“\5 20.0 22.5 

DIVERGENCE AMPLITUDE (l/SEC) *10 8 



LEGEND 

-■ — FEM— A 50 LAT ERS 
- = FEM-A 5 LAYERS 



Figure 66. As in Fig. 62 but for divergence amplitude. 
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c = FDM-A 12 LAYERS 



Figure 67.. As in Fig. 63 but for divergence amplitude. 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR HEATING CASE (FEM-A) 




LEGEND 

= = "EM-A 60 L4TERS 
c = FEM-A 12 LAYERS 



Figure 68. As in Fig. 64 but for divergence amplitude. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR HEATING CASE (FDM-B) 




LEGEND 

= = FDM-B 60 LAYERS 
c = FDM-B 6 LAYERS 



Figure 69. 



As in Fig. 61 but for model FDM-B. 
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HEIGHT. Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR HEATING CASE (FEM-B) 




LEGEND 

a = FEM-B 60 LAYERS 
c = FEM-B 6 LAYERS 



Figure 70. As in Fig. 61 but for model FEM-B. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR HEATING CASE (FDM-C) 




LEGEND 

c = FDM-C 6 LAYERS 
c = FDM-C 60 LAYERS 



Figure 71. As in Fig. 61 but for model FDM-C. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR HEATING CASE (FEM-C) 




TEMPERATURE AMPLITUDE (KELVIN) 



LEGEND 

c = FEM-C 60 LAYERS 
'■ = FEM-C 6 LAYERS 



Figure 72. As in Fig. 61 but for model FEM-C. 



117 



HEIGHT. Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE TOR HEATING CASE (FDM-B) 




TEMPERATURE AMPLITUDE (KELVIN) 



LEGEND 

c = FDM-3 60 LAYERS 
c = FDM-3 12 LATERS 



Figure 73. As in Fig. 63 but for model FDM-B 
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TEMPERATURE AMPLITUDE FOR HEATING CASE (FEM-B) 




LEGEND 

c = FEM-3 60 LAYERS 
c = FEM-B 12 LAYERS 



Figure 74. As in Fig. 63 but for model FEM-B. 
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HEIGHT, Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR HEATING CASE (FDM-C) 




LEGEND 

c = r DM-C 12 LAYERS 
s = FDM— C 60 LATERS 



Figure 75. As in Fig. 63 but for model FDM-C. 
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HEIGHT. Z = - LN(P/PO) 



TEMPERATURE AMPLITUDE FOR HEATING CASE (FEM-C) 




LEGEND 

c = FEM-C 60 LAYERS 
c = FEM-C 12 LAYERS 



Figure 76. As in Fig. 63 but for model FEM-C. 
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boundaries of the heated layer in the consensus profile is 
not reproduced by the 12-layer models, otherwise the profiles 
agree with the consensus outside the heated layer. 

The diabatic heating causes different divergence 
amplitude responses for the B and C grids. The grid C models 
have similar six-level profiles (Figures 77-78) . The minimum 
divergence, which occurs at the center of the heated layer, is 
nearly the same in the six and sixty-layer models. There are 
differences between the two grid C models away from the diver- 
gence minimum. Outside the heated layer the six-layer 
divergence amplitude of the FEM-C model provides the closer 
approximation of the 60-layer consensus profile. Similarly, 
the 12-layer FEM-C profile (Figure 79) is a better approxima- 
tion of the consensus than the 12-layer FDM-C model (Figure 
80) . 

For both six and twelve layers, the grid B models very 
poorly approximate the 60-layer divergence amplitude within 
the heated layer. The 12-layer profiles (Figures 81-82) 
are identical for the grid B models and are closer approxima- 
tions of the consensus 60-layer profile than the six-layer 
results. The 12-layer model is asymmetric within the heated 
layer and does not have the dramatic decrease in divergence 
that exists in the consensus. The twelve and sixty-layer 
profiles are similar outside the heated layer. 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR HEATING CASE (FDM-C) 




LEGEND 

r- = FDM-C 6 LAYERS 
o = FDM-C 60 LAYERS 



Figure 77. 



As in Fig. 61 but divergence amplitude for 
model FDM-C. 
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HEIGHT. Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR HEATING CASE (FEM-C) 




LEGEND 

a = FEM-C 6 LAYERS 
o = FEM-C 60 LAYERS 



Figure 78. 



As in Fig. 61 but divergence amplitude for 
model FEM-C. 
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HEIGHT. Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR HEATING CASE (FEM-C) 
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Figure 79. As in Fig. 63 but divergence amplitude 
for model FEM-C. 
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DIVERGENCE AMPLITUDE FOR HEATING CASE (FDM-C) 




LEGEND 

c = FDM-C 12 LAYERS 
o = FDM-C 60 LAYERS 



Figure 80. As in Fig. 63 but divergence amplitude 
for model FDM-C. 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR HEATING CASE (FDM-B) 




LEGEND 

□ = FDM-B 60 LAYERS 
e = FDM-B 12 LAYERS 



Figure 81. 



As in Fig. 63 but divergence amplitude for 
model FDM-B. 
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HEIGHT, Z = - LN(P/PO) 



DIVERGENCE AMPLITUDE FOR HEATING CASE (FEM-B) 




LEGEND 

" = FEM— 9 60 LATER5 
£ = FEM — B 12 LATERS 



Figure 82. 



As in Fig. 63 but divergence amplitude for 
model FEM-B. 
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IV. CONCLUSIONS 



Although the analytic solutions of each experiment have 
not been derived, some conclusions about the characteristics 
of each model can be drawn. During the Rossby wave and moun- 
tain topography experiments the staggered finite element 
models, FEM-A and FEM-B, display unusual temperature ampli- 
tude behavior in the lowest two layers of the model. The 
oscillation in the temperature profiles of both models may 
be generated by the matrix elements which represent the 
contributions from the eigenfunctions near the surface. It 
is possible that better representations of the staggered 
eigenfunctions near the lower boundary may reduce or eliminate 
this problem. Jagged temperature profiles were not observed 
in the diabatic heating experiment. The 12-hour forecast 
period in the heating experiment may be too short for the 
profile discontinuities to grow to substantial amplitudes. 

The differences between the grids is most apparent in the 
diabatic heating experiment. The 60-layer temperature pro- 
files of all six models represent the same physical solution. 
The FEM-A model has a slightly different 60-layer divergence 
response outside the heated layer than the other five models. 
The grid B and C models have identical temperature amplitude 
profiles for both six and twelve layers. However, the diver- 
gence amplitude profiles are quite different and the grid C 
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models provide the better approximation for both resolutions. 
The grid A models barely respond to the diabatic heating with 
six layers, and poorly approximate their consensus profiles 
with twelve 'layers. There are indications that both grid A 
models may not converge to their 60-layer temperature and 
divergence solutions with much less than 60-layer resolution. 

In all experiments, the grid B model produces virtually 
identical results, aside from the small amplitude oscillation 
in the FEM-B temperature profiles. Based on these experiments, 
no accuracy is gained by using the finite element approxima- 
tions with grid B. In two of the three experiments the grid 
C results are identical. The FEM-C model provides a closer 
approximation of the consensus 60-layer divergence profile 
than the FDM-C model in the diabatic heating experiment. 
Generally, the lower resolution FEM-A model does not approxi- 
mate its 60-layer solution as well as the FDM-A model. It is 
not known if the different approximation characteristics of 
the FEM-A model are due solely to its profile problems at the 
bottom of the atmosphere. 

The importance of the -^ ( Z ) and eigenfunctions in 

the staggered finite element models is not clear. They are 
prominent terms in the mathematical derivations, but their 
physical significance is not apparent. 
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APPENDIX A 



FINITE DIFFERENCE APPROXIMATIONS 



ciu 

1 . For terms of the form -5-=- W: 

dZ 

a. FDM-A and FDM-B, at level Z = Z ± : 



du 

dz' 



i u. , . -u. u. - u. , 

« - K-'r 1 •' 1_1 ’ 



i+l 



— u . u . — U • _ 1 

- Z . } + w i-l* ( Z . - Z . , ^ 



i-1 



b. FDM-C, at level Z = Zj: 



du 



dZ ’ W 2 ' W i‘ [ ( Z ! 



u, . - u . 

+ ( 

i+l Z i 



u . - u . . 
1 l-l 

Z ! - Z ! , 
1 l-l 



) ] 



3 T 

2 . For terms of the form W: 

aZ 

a. FDM-A, at level Z = Z^: 



3 T 
3 Z ‘ 



W = 



2 [ V<Z. 



T. -t. 
1+1 1 



i + l - z i 



) + W. , • ( 



i-i > Zi -z^' 



b. FDM-B and FDM-C, at level Z = Z.: 

1 



= — *W • [ ( - .f + - L 1 ~ 1 ) + ( 1 1 i.) l 

3Z 2 i U Z ' - Z ! ' ! - Z ! * J 



T - T 

i+l i-1, 



T. -T. 



i-1, 



i+l 



i-1 



du 

3 . For terms of the form 3-=- V: 

dZ 

a. FDM-A, at level Z = Z-: 

1 



du 

dZ* 



u . - u . u.-u. , 

V = i-V[(— -) + (— — 

2 L V Z . . - Z . ; K Z . - Z . . 

l+l 1 1 l-l 



) ] 
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b. FDM-B, at level Z 



Z i 



du 
dZ ’ 



V 




+ V i ) '(^ ±i 
1 Z i+1 




c. FDM-C, at level Z = Z[: 



du 

dZ 



1 V 



|*v* [ ( "f +1 gt ) 

2 z i+l _z i 



u . - u . , 
(__i X -~) 1 

l Z ! - Z ! , ' J 

• l l-l 
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APPENDIX B 



GALERKIN FORM OF FEM-C PROGNOSTIC EQUATIONS 



1. Vorticity Equations (2-35) - (2-36) : 



i+1 dA 



3 Zr 



l '^zr- / 4>i (Z)4>i (Z)dZ = 

i=i-l ar z„ J 



i+1 , T 

-f l D? (t) / <0 . ( z) 4). (z)dz 

j=i-i z 0 3 



i+1 i+1 



- .k 



V l l u.A,(t) / 4> . (Z) <L (Z)4>. (Z)dZ 

j = i-l k=i-l 11 * z n 0 K 1 



i+1 • T 

6 l vut) J <p.(z)<p.(z)dz 
j=i-l z„ 3 



(B-l) 



i+1 dA^ Z T 

l '/ ^ (Z) 4-i (Z)dZ = 

j=i-l Z 0 ] 



i+1 



-f l d 2 (t) / (j) . (Z) 4). (Z)dZ 

j=i-l z n 3 



i+1 i+1 



+ V l 

j 



l l u.A^(t) / 4> . (Z) 4> k (Z) 4). (Z) dz 

i=i-l k=i-l - 1 Z„ J 



i+1 • T 

6 l v2(t) / 4> . (Z) 4>. (Z)dZ 

j=i-l z 0 3 



( B-2 ) 
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a. Divergence Equations (2-37) - (2-38) : 



i +1 

I 

j=i-l 



i +1 

l 

j=i"l 



dD? Z T i+1 • Z T 

-rr± / (J). (Z) 4>. (Z)dZ = f l A?(t) / <$> ■ ( Z ) 4 > . ( Z ) dZ 

z o ’ j' 1 ' 1 z 0 



i +1 i +1 



• “ V 



- _k 



l l u .D« (t) / 4), (Z)cJ>, (Z)<f>. (Z)dZ 

j=i-l k=i-l 3 Z Z. 3 K 1 



i+1 • T 

B l U? (t) / <$>. (Z) <|> . (Z)dZ 

j=i-l z Q ] 



i +1 i +1 



- „k 



'T d 4> . 



y l l u,w (t) / <}). (z) 41. (z) dz 
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_ i+i , ' T 

+ y l H ; 1 (t) / <J>. (Z) <j> 

i=i-l x Z rt 3 1 



(Z) dZ 



(B-3) 
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Z 0 



/ <j). (Z) 4> i (z)dz = f _ l A^(t) / <|). (z) 4> i (z)dz 



j=i-l 



i +1 i +1 , *T 

+ V l l u.D*(t) / 4>. (Z) 4> (Z) 4 ,. (Z)dZ 

j=i-l k=i-l 3 1 Z : K 1 
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- B I Uo(t) / 4 >. (Z) 4 >. (Z)dZ 

i=i-l Z„ J 



_ k - J<r d( ^ 
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+ v l l n u i w 1 (t) / ^ 4> k (z)4> i (z)dz 

■o 
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i+1 . T 

^2 (t) / 

j=i-l z Z, 



+ M 2 _ J , H 2 (t) / <f>j (Z) 4>i(Z)dZ 



(B-4 ) 
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3. Thermodynamic Equations (2-39 )- (2-40 ) : 



i+1 dT^ Z T 



l j7~ / 4 >. (z) 4 > . (z) dz 

i=i-l J 



i+1 i+1 



y l l u.T^(t) / <t > . (Z)cD , (Z)<j> (Z)dZ 
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f i+1 i+ l _ k T d<j>. 

+ <f> n u j v i (t) 7 / dz 2 *k (z) *i (z)dz 

3 = 1-1 k=i-l J Z Q 



i+1 i + 1 , T d(J> . 

.1.1 T j w i (t) / ar 1 *k (z) n (z)dz 

]=i-l k=i-l J Z Q 



i+1 j T 



+ l Q^(t) / 4>, (Z)<|» i (Z)dZ 
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i+1 dT^ Z T 
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j-i-l Zn 3 



i+1 i+1 



- Jc 
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+ y Q, (t) / <J> • ( Z ) <J) . (Z) dZ 
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APPENDIX C 



MATRIX ELEMENTS FOR FEM-C 



1. Notation: A^ = Z^ but A 2 = Z 2 ~ Z i 



A! = Z! -Z! , 
i 1 l-l 



u . , T . 

l l 



u^ T at i 



th 



vertical level 



2. Matrix definition: M (row, column) 

th 

The i row of a matrix contains the weighting factors of 

the terms that affect the value of the forecast variable at 

the i th vertical level. The weighting factor in the j 1 "* 1 

th 

column of the i row is the influence on the forecast varia- 
ble at the i^ level from the variable at the j ^ vertical 
level. Level i = 1 is the surface, level i = N+l is the top 
of the atmosphere. 

3. To calculate the terms of each matrix, substitute in the 
matrix formula for 4>^ + ^(Z), <J>^(Z), 4>^_^(Z) using functions 
defined in local coordinates with respect to <j>^(Z). Define 
S = Z -Z[. 

The eigenfunctions in local coordinates are 



CO = 



E, + A! 

l 

A! 

l 

- 5+A bi 



-A| <_ C £ 0 



0 1 ^ 1 A i + 1 



(C-l ) 
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4 > 



i+1 



(€) 




0 1 ^ i A! + i 



(C-2 ) 






- r 




-A! < £ < 0 . 



(03) 



To calculate the matrix elements corresponding to the surface 
and the top of the atmosphere, use the above formulas for 
i = 1 and i = N+l, respectively. 

4. Mass matrix, M: 

a. General formulas, i = 2,...,N: 

M ( i , i-1) = |A i (04) 

M ( i , i ) = |(A i + 1 + A i ) (05) 

M ( i , i+1 ) = |a 1 + 1 (06) 

b. At level i = 1: 

M ( 1 , 1 ) = -|a 2 (07) 
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M ( 1 / 2 ) 



(C-8 ) 



M ( 1 / 2 ) = 


ifi 2 (C-8) 


c. At level i = 


N+l: 


M(N+1,N) = 


■ ! a n + i (c - 9 » 


M (N+l , N+l) 

V 

5. Matrix N: 


= sVi (c - l0 » 



a. General formulas, i = 2,...,N 



N ( i , i-1) = 


= 1 ^A i - (u i +u i _ 1 ) (Oil) 


N ( i , i ) 


¥“i' u i+i + A i> + n" ( “i + i' A i + i 



+ u i _ 1 -A i ) (012) 



N ( i , i+1) = 


12 A i+l' (u i+l + u i ) (C 13) 


b. At level i = 


1: 


N ( 1 , 1 ) = 


T U 1 * A 2 + 12 U 2* A 2 (014 ) 


N ( 1 , 2 ) = 


Yj- A 2 - (u 2 +u 1 ) (015) 


c. At level i = 


N+l: 


N (i , i-1) = 


= ^-Ai ( Ui+Ui_i) (016) 


N(i,i) = 


rv^i + ir V“i-i (c " 17) 
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6. Matrix P: 

a. General Formulas, i = 2,...,N 



P (i, i-1) = 



l (u i - u i-r 



(C-18) 



P(i,i) = 



| (u 1+ l -“i-i) 



(C-19 ) 



P (i, i+1) = 



l (u i+i - u i> 



(C— 20) 



b. At level i = 1: 



P(l,l) = 



I (U 2 "V 



(C-21) 



P ( 1 / 2 ) = 



| (U 2 - U l } 



(C-22) 



c. At level i = N+l: 



P (N+l , N) 



- u. T ) 
6 N+l N 



(C-23 ) 



P (N+l , N+l ) 



I (U N + 1 ‘V ' 



(C-24) 



7. Matrix R: 

a. General formulas, i = 2,...,N 



R(i,i-1) 



?< T i 



(C-2 5 ) 



R(i,i) = 






(C-26) 



R(i , i+1) = 



i (T i + i - T i> 



(C-2 7 ) 
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b. At level i 



1 



R ( 1 / 1 ) = -|(T 2 -T 1 ) (C-28 ) 

R(l,2) = |(T 2 -T x ) ( C— 2 9 ) 

c. At level i = N+l: 

' R (N+l , N) = ^(T n+] _ -T n ) (C-30) 

R (N+l , N+l ) = |(T n+] _ -T n ) (C-31) 
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APPENDIX D 



GALERKIN FORM OF FEM-A PROGNOSTIC EQUATIONS 

1. The vorticity equations, (2-35) - (2-36) , have the same 
form as the vorticity equations in FEM-C. While the FEM-A 
equations look the same as equations (B-l)-(B-2) in Appendix 

V 

B, the equations for <j), , <j) . and <j). are different because the 

1 J K 

<f>(Z) eigenfunctions are defined for the staggered levels in 
model FEM-A. 



2. Divergence Equations (2-37) - (2-38) : 



i+1 dD^ Z T i+1 Z T 

l ^ / 4). (Z) 4). (Z)dZ-= f J A? ( t ) / 4) . (Z) 4> • (Z)dZ 

:=i-i z 0 3 j=i-l z 0 3 

i+1 i+1 _ k Z t 

- y [ l u.D (t) / 4>. (Z) 4, (Z) <p. (Z)dz 

j=i-l k=i-l 3 z Z 3 K 1 



i+1 • T 

-6 l u 3 (t) / <p . (Z)<p. (z)dz 
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i+1 i+1 



-y l l u.w£(t) / ^4 4* k (z)4,.(z)dz 

j = i-l k=i-l 3 Z Z. ciz K l 



J T d4> . 



2 i+1 i T 

+ y J (t) / 4, . (z) 4). (Z) dz 

j=i-l 1 Z 0 3 1 



(D-l) 
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i+1 dD^ Z T 



i+1 



I jl / 4> i ( Z ) (Jy- ( Z ) dZ = f l A~ ( t ) / <t> . (Z) <j>. (Z)dZ 
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(D-2) 



3. Thermodynamic Equations (2-39 )- (2-40) : 



i+1 dT 3 Z T 



I / 4> i (Z)<j>.(Z)dz 

j=i-l dt z n 3 
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= - y l 
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i z 

i+1 dT^ T 



l 3T— / <z> 4> ± (Z)d2 

j=i-l dt z. ] 



i+1 i+1 



- „k 



y - I I u.T (t) / <|).(z)<o k (z ) 4 >. (z) dz 

i=i-l k=i-l 31 Z rt 3 



i+1 i+1 



T dcf> . 



+ ^ 



<I> . I „ l ,vS (t) , / ar 1 * k (z) i’i (z)dz 

„ 3=1-1 k=i-l J Z 



i+1 i+1 



T d<j> . 



J.TJ. J.TJ. _ 1 , ± '-I'f ’ 

. j , J j w 2 (t) / ar 1 (z) +i (z) dz 



i+1 



+ I Qn ( t) / 

j=i-l z n 3 



(Z) dZ 
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APPENDIX E 



LOCAL COORDINATES FOR STAGGERED <t> ± ( Z ) 

1. To evaluate the integrals in the finite element approxima- 
tions of the vorticity and divergence equations in models 
FEM-A an£ FEM-B, express <j>^ / ^i+i' ^i' ^i-1 

and ib . 0 in terms of the local coordinate E = Z - Z . . The 

r i-2 ^ 1 

general form of the local coordinate system is shown below. 




2. In terms of local coordinates for 4 k(£), the equations 

for <t . , , A . , <p . . , v . i » A • > ip ■ i and ib . _ are listed below . 
Y i+1 l *i-l l+l Y i l-l +-2 

The notation used is defined in Appendix C. 
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